Ab initio Monte Carlo simulations for finite-temperature properties: Application to 

lithium clusters and bulk liquid lithium 
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Ab initio Monte Carlo simulations have been performed to determine the equilibrium properties 
of liquid lithium and lithium clusters at different temperatures. First-principles density-functional 
methods were employed to calculate the potential-energy change for each proposed change of config- 
uration, which was then accepted or rejected according to the Metropolis Monte Carlo scheme. The 
resulting structural properties are compared to data from experimental measurements and ab initio 
molecular dynamics simulations. It is shown that accurate structural information can be obtained 
with ab initio Monte Carlo simulations at computational costs comparable to ab initio molecular 
dynamics methods. We demonstrate that ab initio Monte Carlo simulations for the properties of 
fairly large condensed-matter systems at nonzero temperatures are feasible. 
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I. INTRODUCTION 

Molecular dynamics (MD) and Monte Carlo (MC) 
simulations have been the major techniques for calcu- 
lating finite-temperature properties of condensed-matter 
systems. In particular, ab initio MD simulations, 
which combine classical molecular dynamics with first- 
principles quantum theory, determine the atomic and 
electronic structures of a system simultaneously and ac- 
curately, as well as the temperature-dependent proper- 
ties and time evolution of the systemiiiS Ab initio MD 
methods have been widely used to provide reliable and 
accurate information for various systems since Car and 
Parrinello first introduced their work in 1985ii 

On the other hand, ab initio Monte Carlo simulations, 
which employ first-principles quantum theory to calcu- 
late the potential energy of a system at each classical 
MC step, have so far been used only to a very limited 
extent. The majority of MC simulations for condensed- 
matter systems have been performed with empirical or 
semiempirical atom-atom interaction potentials, which 
are fitted to available experimental data or ab initio cal- 
culations. While the atom- atom potentials may work 
well in some cases, they do not provide reliable results 
in other instances. Since the quantum many-body ef- 
fects and the detailed electronic structures cannot be 
included in empirical or semi-empirical potentials, first- 
principles theory would obviously be desirable for provid- 
ing reliable energetics. The major concern with ab ini- 
tio Monte Carlo simulations is the high computational 
demand. In particular, Monte Carlo simulations with 
empirical potentials often generate millions of sampling 
configurations. It is a formidable task to determine the 
potential energies of such a huge number of configura- 
tions with ab initio total-energy calculations. In recent 



years, however, several groups have performed ab ini- 
tio MC simulations for very small clusters, and they 
have tried to reduce the number of sampling configura- 
tions (to typically io4), 3,4,5,6,7,8,9,10,11, 12,13,14 ighjkawa et 
^j3ji2ji4 appijg(;j g^f) ifiHio calculations to the Monte Carlo 
simulated annealing algorithni^^. They performed sim- 
ulations for the Lig cluster, small hydrogenated lithium 
clusters (LisH, Li4H2, and LijH), and HCl(Il20)„ (n = 
3, 4) clustersAi2ii4 Other groups"'-'^'^'-'^" used the ab ini- 
tio molecular orbital (MO) method in classical MC sim- 
ulations to study the interactions of ions and molecules 
with small water clusters. Structural and thermal prop- 
erties of water dimersii and lithium clusters (Li4, Li^, 
and Li8)Si24i^ were also investigated with ab initio MC 
simulations. 

We emphasize that all the previous ab initio MC sim- 
ulations were limited to very small clusters. Since exper- 
imental data for the structures of the clusters at nonzero 
temperatures were not available, direct comparisons of 
the above simulations with experiments were not possi- 
ble. Therefore, the accuracy and efficiency of ab initio 
MC simulations remain to be determined. Furthermore, 
the feasibility of the extension of ab initio MC simula- 
tions to bulk systems, which require models containing 
many more atoms than the small clusters, has yet to be 
tested. As a first application of ab initio MC simulations 
to bulk systems, we have chosen the liquid phase of bulk 
lithium. We have also investigated the structural proper- 
ties of a cluster containing 16 lithium atoms at different 
temperatures. The simulated results for the lithium clus- 
ter are physically reasonable, and the structure of liquid 
lithium obtained from our ab initio MC simulations is 
in good agreement with the available experimental mea- 
surements. We also demonstrate that the computational 
demands of simulating liquid lithium with our ab initio 
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MC algorithms are comparable to those of ab initio MD 
simulations. 

The remainder of this paper is organized as follows. 
Section II describes our ab initio MC algorithm. Section 
III describes the simulation results for a 16-atom lithium 
cluster (Sec. Ill A) and for liquid bulk lithium (Sec. Ill 
B), with a discussion of general aspects and possible im- 
provements of our ab initio MC algorithm in Sec. Ill C. 
Sec. IV contains our conclusion. 



II. AB INITIO MONTE CARLO ALGORITHM 

To simulate the cluster of 16 lithium atoms, the atoms 
were put in a periodically repeated cubic supercell. The 
size of the supercell (10 x 10 x 10 A'^) was much larger 
than the concentrated cluster (see Fig. 1.). Bulk lithium 
was represented by a cubic box (supercell) with peri- 
odic boundary conditions containing 54 Li atoms. The 
size of the supercell (10.698 x 10.698 x 10.698 A^ and 
10.874 X 10.874 x 10.874 A^ for simulations at T = 523 K 
and T — 725 K, respectively) was determined to match 
the experimental density of liquid lithium at the given 
temperaturepi^ 

The simulations were performed in canonical ensem- 
bles. The standard Metropolis algorithm, ^^ -'^^ which pro- 
vides an efficient approach for simulating the equilibrium 
properties of an atomic system at a given temperature, 
was used in our ab initio MC simulations. Starting from 
a pre-selected initial configuration, the simulations re- 
peated the following steps. (1) Given the current con- 
figuration / of the system, a new configuration J is gen- 
erated by making a small random displacement of one 
randomly chosen particle (suitable for small systems) or 
random displacements of all the particles simultaneously 
(more efficient for systems containing a large number of 
particles). (2) Once the total energies {Ej and Ej) of 
these two configurations are calculated, the acceptance 
probability of the new configuration J is then determined 



P{J\I) = min 



1, exp 



-jEj-Ei) 



(1) 



where T is the temperature and ks is Boltzmann's con- 
stant. (3) If the configuration J is accepted, it serves 
as the current configuration of the next MC step, and 
El is set equal to Ej. If the configuration J is not ac- 
cepted, the configuration / and its energy Ej are retained 
and used to start the next step. In this way, the sys- 
tem will eventually reach equilibrium and evolve toward 
a Boltzmann distribution>iii^ The resulting sequence of 
sampling configurations is then used to obtain the equi- 
librium properties of the system. 

To calculate the total energy at each MC step, we em- 
ployed an ab initio pseudopotential density-functional- 
theory (DFT) method with a plane- wave basis set. 
All the ab initio total-energy calculations were car- 
ried out using the Vienna ab-initio simulation package 



(VASP)ii2422i2i The electronic exchange-correlation ef- 
fects were treated with the generalized gradient-corrected 
functional given by Perdew and Wang i^^'^? We adopted 
the Vanderbilt ultrasoft pseudopotential to replace the 
core electrons^^-^'"' and used the conjugate-gradient tech- 
nique for performing electronic relaxations^. The total- 
energy calculations were conducted with 4 k-points in 
the three-dimensional Brillouin zone with a plane-wave 
energy cutoff of 120 eV. The convergence of the total 
energy differences between different configurations was 
checked for selected configurations with higher cutoff en- 
ergies (up to 200 eV) and more k-points (up to 16 k- 
points), and the convergence was found to be within a 
few meV. Because we used a pseudopotential approach, 
the total-energy calculations for the configurations with 
too small distances between atoms, which result in core- 
core overlap, would give wrong total energies. We there- 
fore defined a hard-core distance of 2.0 A between any 
two Li atoms. Trial configurations with any interatomic 
distance smaller than the cutoff distance were always re- 
jected in the MC simulations. This is physically reason- 
able because such configurations would have very high 
potential energies due to the strong core-core repulsion. 



III. RESULTS AND DISCUSSION 
A. Small cluster 

Starting from the configuration with the atoms at or- 
dered bcc positions (shown in Fig. 1), wc first performed 
our ab initio MC simulations for the cluster of 16 lithium 
atoms at temperatures of 1000 K and 5000 K. The MC 
moves in the simulations for the cluster consisted of ran- 
dom displacement of one randomly selected lithium atom 
at each MC step, with the maximum displacement ad- 
justed to give an overall acceptance rate of about 50%. 
The chosen atom was displaced from its old position with 
equal probability to any position inside a sphere sur- 
rounding the old position, with the radius of the sphere 
being the maximum displacement. We obtained maxi- 
mum displacements of 0.75 A and 1.2 A for the simula- 
tions at the temperatures of 1000 K and 5000 K, respec- 
tively. About 10000 MC steps were needed to reach the 
equilibrium state. Averages for calculating the structural 
properties were taken over the next 10000 sampling con- 
figurations after equilibration. It took approximately 1 
minute of CPU time to obtain the total energy of each 
configuration with a single 375 MHz processor on an IBM 
SP3 supercomputer. 

The pair correlation function g{r) can be defined asS^ 



2V 



1 



N{N - 1) 47rr2 



(EE^(^-^'^))' (2) 



'j<i 



where V is the system volume, N is the number of parti- 
cles, Tij is the distance between a pair of particles i and j, 
and (•) signifies averaging over the sequence of simulated 
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FIG. 1: The initial configuration of the cluster of 16 lithium 

atoms for the ab initio MC simulations. The distance be- 
tween two nearest-neighbor atoms is 2.96 A, the same as that 
determined by our ab initio total-energy calculations for the 
crystalline bulk Li. The size of the supercell is 10 x 10 x 10 

A3. 



configurations. The pair correlation function is propor- 
tional to the probability of finding a pair of particles a 
distance r apart, normalized by the square of the parti- 
cle density such that it approaches unity in the limit of 
a uniform random distribution at the same density. It 
hence provides a measure of the local spatial ordering of 
the system. As described in Ref. 26, wc determined g{r) 
by dividing the range of r ([2 A , 10 A]) into 400 intervals 
of length 0.02 A (Ar), calculating every Vij for a given 
configuration, and recording in narrow bins [r, r + Ar) 
the frequency of occurrence of the different particle sepa- 
rations. In this way a histogram was built up and finally 
normalized at the end of the scan over the configurations. 

Figure 2 shows the pair correlation functions for the 
cluster at temperatures of 1000 K and 5000 K, respec- 
tively. It can be observed that there are two peaks in the 
pair correlation function at T = 1000 K. This indicates 
the existence of two atomic shells for the structure, and 
therefore the system locally has a liquid-like structure. 
At a temperature of 5000 K, however, there is no obvi- 
ous second peak in the pair correlation function, and thus 
the system has a gas-like structure. Given that the boil- 
ing point of bulk Li is about 1600 K, the results from our 
simulations for the cluster of lithium atoms are physically 
reasonable. Note that the pair correlation function g{r) 
at T = 5000 K does not approach zero when r is close 
to 2.0 A, indicating that the configurations with very 
small separations (smaller than 2.0 A) between atoms 
would contribute to the ensemble averages. This reflects 
the gas-like structure of the system at T = 5000 K be- 
cause the probability of occurrence of configurations with 
very small separations between particles is expected to 
be higher for a gas-like structure than for a liquid-like 




2.0 4.0 6.0 8.0 10.0 



r(A) 

FIG. 2: Equilibrium pair correlation functions g{r) for a clus- 
ter of 16 lithium atoms at temperatures of 1000 K and 5000 K, 
obtained from our ab initio Monte Carlo simulations. 



structure. 



B. Bulk system 

Following the simulations for the lithium cluster, we 
carried out ab initio MC simulations for the liquid phase 
of bulk lithium. In contrast to the simulations for the 
lithium cluster, the Monte Carlo moves in the simulations 
for liquid lithium consisted of random displacements of 
all the atoms in the supercell simultaneously. The other 
aspects of the MC moves were the same as in the sim- 
ulations for the lithium cluster. In order to achieve an 
overall acceptance rate; of roughly 50%, we found that the 
maximum displacements at temperatures of 523 K and 
725 K were 0.05 A and 0.06 A, respectively. We started 
our simulations at T = 725 K from an initial configu- 
ration with all the Li atoms at their bulk bcc positions. 
Roughly 8000 MC steps were required to reach equilib- 
rium. Since the initial configuration was a perfect crys- 
talline structure, far away from the equilibrium structure 
of the liquid phase of bulk lithium, it required relatively 
long MC runs to reach equilibrium. In contrast, the ini- 
tial configuration for the simulations at T = 523 K was 
chosen as one of the equilibrium sampling configurations 
obtained in the simulations at T = 725 K. A reduced 
number of MC steps (5000 MC steps) was then required 
for the equilibration period. After equilibration at either 
temperature, we generated 5000 sampling configurations 
for averages. 

Our simulations were carried out in parallel on four 
375-MHz processors of an IBM SP3 supercomputer. The 
total-energy calculation at each MC step required ap- 
proximately 150 seconds of CPU time. On average, 
the electronic relaxation converged within seven self- 
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FIG. 3: Equilibrium pair pair correlation function g{r) for 
liquid lithium at a temperature of 523 K, calculated from ab 
initio Monte Carlo and ab initio molecular dynamics simula- 
tions (histograms). Also shown are experimental data mea- 
sured by X-ray diffraction (circles) 



consistency steps. The software (VASP) that we used 
for the total-energy calculations contains computations 
for the forces after the electronic relaxation converges. 
Such computations are not needed in ab initio MC sim- 
ulations. Should VASP be modified to exclude the force 
calculations, the computational costs would therefore be 
decreased by approximatively 10% - 15%. 

For comparison, we also performed ab initio MD 
simulationsiSiSfliSi for liquid lithium. Canonical ensem- 
bles were generated with the algorithm of Nose,^'' and the 
time step was set to 1 fs. The initial configuration for the 
ab initio MD simulations at both T = 725 K and T — 
523 K was the same perfect bcc structure as that for the 
ab initio MC simulations at T = 725 K. The equilibration 
periods extended over 5000 and 3000 time steps at T = 
523 K and T = 725 K, respectively. Ensemble averages 
were calculated from the sampling configurations of the 
next 5000 steps. Each step required about 400 seconds 
of CPU time on an IBM SP3 supercomputer (in parallel 
on 4 processors), and the corresponding wall-clock time 
was roughly 100 seconds. 

The pair correlation functions for liquid lithium at 
temperatures of 523 K and 725 K, obtained from our 
ab initio MC and ab initio MD simulations are shown 
in Fig. 3 and Fig. 4, respectively. For comparison, the 
experimental data for liquid lithium at T = 523 K ob- 
tained by X-ray diffraction^S, are also shown in Fig. 3. 
We observe that the agreement between ab initio MC, 
ab initio MD, and experiment is very good. In particu- 
lar, the three peaks determined from our ab initio MC 
simulations are essentially located at the experimental 
positions. The intensities of the peaks obtained from 
the ab initio MC simulations are also basically in agree- 



FIG. 4: Equilibrium pair correlation function g{r) for the 
liquid phase of lithium at a temperature of 725 K, obtained 
from ab initio MC and ab initio MD simulations. 

ment with those from both the ab initio MD simulations 
and the X-ray data. The small differences of the inten- 
sities of the peaks between the theoretical results and 
the experimental data may be due to the small size of 
our supercells. The experimentally observed tempera- 
ture effects,-2S that is, broadening of the peaks and slight 
shifts to higher positions of the second and third peaks at 
the higher temperature are also well reproduced by our 
ab initio MC simulations (see Fig. 3 and Fig. 4). 



C. General discussion 

Finally, we discuss some general aspects and possible 
improvements of our ab initio MC algorithm. The com- 
putational costs of ab initio MC simulations are domi- 
nated by the ab initio total-energy calculations. Since 
efficient iterative schemes for ab initio total-energy cal- 
culations with a plane-wave basis set have achieved order 
iV^ (iV is the number of the atoms in the supercell) scal- 
ing for systems containing up to 1000 electrons?^ the cost 
of simulations with our ab initio MC algorithm would in- 
crease essentially as the square of the number of atoms. 

While the computational costs are comparable, our ab 
initio MC algorithm is not any faster than ab initio MD 
simulations. In particular, a larger number of steps are 
required to reach equilibrium in our ab initio MC sim- 
ulations for bulk lithium than in the corresponding ab 
initio MD simulations. For the lithium systems that we 
have investigated, our ab initio MC algorithm produces 
the required sampling configurations at a speed similar 
to the ab initio MD simulations. In this respect, our ab 
initio MC simulations are as efficient as ab initio MD 
simulations. 

In some cases, like simulations for strongly covalent 
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systems and situations involving volume changes, how- 
ever, our ab initio MC algorithm may not be efficient 
due to the low acceptance rate of the traditional MC 
moves. In such situations, therefore, our ab initio MC 
algorithm needs to be improved. A cluster MC move 
scheme proposed by Lee and Swendsen^ii may be used 
for this purpose. We have tested the cluster MC move 
scheme for the bulk lithium systems that we investigated 
with traditional MC moves as described in Sec. III. We 
found that the use of cluster MC moves in ab initio MC 
simulations provided almost the same performance as the 
traditional MC moves, indicating that both kinds of MC 
moves might work equally well for metallic systems. It 
is also possible that the relative performance of different 
kinds of MC moves somewhat depend on the number of 
particles included in the supercell. It would be interest- 
ing to see if the cluster MC move scheme would be more 
efficient than the traditional MC approach in simulations 
for systems containing large numbers (for example, hun- 
dreds) of particles. In the cases where the boundary con- 
ditions or/and the volume of the simulated cells change, 
which result in a low acceptance rate with the traditional 
MC moves, the cluster MC move scheme has been shown 
to increase the acceptance rate significantlyi^ We also 
expect that the cluster MC move scheme would be more 
efficient than the traditional MC move method in sim- 
ulations for systems containing strong covalcnt bonds. 
Another approach, a hybrid of MD and MC simulation 
in which trial MC moves are generated with classical MD 
methods, may also enhance the efficiency of ab initio MC 
simulations^ Other techniques aiming at sampling im- 
portant configurations in the configuration space include 



the J-walking procedure proposed by Frantz et ali**i2iS 

and a big-move method suggested by Akhmatskaya et 
alA 



IV. CONCLUSION 

In conclusion, we have shown that the use of ah initio 
calculations in classical Monte Carlo simulations provides 
structures of liquid lithium in agreement with experimen- 
tal measurements and ab initio MD simulations. Our 
work also demonstrates that ab initio MC simulations 
for condensed-matter systems with up to several tens of 
atoms (and possibly with more than a hundred atoms) 
are feasible and may provide an alternative way to inves- 
tigate finite-temperature properties of such systems. 
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